function [celltab,xsim]=plotgamma(alpha,beta,into); 
%function [xsim,celltit]=plotgamma(alpha,beta,into); 
% Plots the pdf of an B(alpha,beta) distribution 
%   including mean and standard deviation
% Legend provides percentiles 
% into.min    minimum value to plot
% into.max    maximum value to plot
% into.step   step between points
% into.save   =1 save 
% into.fpath  file path 
% into.opath  output path
% into.per
% Output
% -------
%  celltit   cell, contains percentiles 
%  xsim      vector of simulated values used to compute percentiles 
%
% Modified 6/02/03
% ------------------------------------------------------------------
meant=alpha*beta; 
vart=alpha*(beta^2);
sd=sqrt(vart);
if nargin < 3  | isempty(into) == 1
    into.min=0.01; 
    into.max=15; 
    into.step=0.05; 
    into.save=0; 
end
x=(into.min):(into.step):(into.max); 
pdf=gampdf(x,alpha,beta); 
figure; 
%subplot(2,1,1); 
plot(x,pdf); 
XLim([into.min-2*(into.step) into.max+2*(into.step)]);
astr=num2str(alpha); 
bstr=num2str(beta); 
mstr=num2str(meant); 
sstr=num2str(sd); 
t1=['Gamma(' astr ',' bstr ')   mean : ' mstr ' sd: ' sstr ]; 
title(t1);
% Simulate 
xsim=gamrnd(alpha,beta,20000,1); 
per=[0.001 0.01 0.05 0.1 0.5 0.9 0.95 0.99 0.999]; 
perc=num2cell(per);
perc2=gaminv(per,alpha,beta); 
perc2=num2cell(perc2'); 
jj=1; tstr=length(perc); 
celltit=cell(tstr,1); 
celltab=cell(tstr,2);
for jj=1:tstr; 
    celltit{jj}=[perc{jj} '    ' perc2{jj}];
    celltab{jj,1}=perc{jj}; 
    celltab{jj,2}=perc2{jj};
end 
legend('density',celltit{1},celltit{2},celltit{3},celltit{4},celltit{5},celltit{6},celltit{7},celltit{8},celltit{9})
celladd=[{'1'} {mstr}; {'2'} {sstr}]; 
celltab=[celladd;celltab]; 
if into.save==1 
    into.desc=['Gamma' astr '_' bstr]; 
    savegraph(into)
end